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Abstract 

Linear response of simple (i.e., condensed) Bose-Einstein condensates is known to lead to the 
Bogoliubov- de Gennes equations. Here, we derive linear response for fragmented Bose-Einstein 
condensates, i.e., for the case where the many-body wave function is not a product of one, but of 
several single-particle states (orbitals). This gives one access to excitation spectra and response 
amplitudes of systems beyond the Gross-Pitaevskii description. Our approach is based on the 
number-conserving variational time-dependent mean field theory [O. E. Alon, A. I. Streltsov, and 
L. S. Cederbaum, Phys. Lett. A 362, 453 (2007)], which describes the time evolution of best- 
mean field states. Correspondingly, we call our linear response theory for fragmented states LR- 
BMF. In the derivation it follows naturally that excitations are orthogonal to the ground-state 
orbitals. As applications excitation spectra of Bose-Einstein condensates in double-well potentials 
are calculated. Both symmetric and asymmetric double-wells are studied for several interaction 
strengths and barrier heights. The cases of condensed and two- fold fragmented ground states are 
compared. Interestingly, even in such situations where the response frequencies of the two cases 
are computed to be close to each other, which is the situation for the excitations well below the 
barrier, striking differences in the density response in momentum space are found. For excitations 
with an energy of the order of the barrier height, both the energies and the density response of 
condensed and fragmented systems arc very different. In fragmented systems there is a class of 
"swapped" excitations where an atom is transfered to the neighboring well. The mechanism of its 
origin is discussed. In asymmetric wells, the response of a fragmented system is purely local (i.e., 
finite in either one or the other well) with different frequencies for the left and right fragments. 
This finding is in stark contrast to that for condensed systems. 

PACS numbers: 03.75.Kk, 05.30.Jp, 03.65.-w 
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I. INTRODUCTION 



The precision with which ultra-cold bosonic quantum gases and Bose-Einstein conden- 
sates (BECs) can be manipulated nowadays [H |2] has allowed to study not only their ground 
states but also excited states and dynamics, such as dynamical splitting of a BEC [3], Joseph- 
son dynamics [U |3] , dynamical creation of number squeezed states [HI E] , quantum optimal 
control [8] , and multi-band physics P, . 

Crucial for the understanding of the dynamics of quantum gases are excitation spectra. 
Theoretically, they have been widely studied for BECs using the standard Bogoliubov-de 
Gennes (BdG) equations, which can be regarded as the linear response of the ground state 
described by the Gross-Pitaevskii (CP) equation [TTHT^ to an external perturbation. Pre- 
dicted spectra [ISHIZj compared accurately to experiments measuring for example collective 
excitations of trapped BECs [TSl [H] , the speed of sound propagation of quasi-particles [20] , 
and excitations of BECs in the bulk regime [211 ESj. Similar theoretical and experimental 
studies have been performed for two-species BECs, modeled by two coupled Gross-Pitaevskii 
equations (2GP) and the corresponding Bogoliubov-de Gennes equations [23H27j. 

In those works it is assumed that the atoms are essentially condensed. Therefore, the 
ground state of the system is well described by the CP (or 2GP) equation. However, there 
are numerous examples where this is not the case and, instead, one deals with fragmented 
condensates [2HHS2]- Examples involve condensates in symmetric and asymmetric 

double- well potentials [301 EH ES] • Mott insulator states in few- well systems [5S1 SD] and op- 
tical lattices [SI US] represent multiple fragmented states. Other examples of fragmentation 
can be found in cold atom systems exhibiting translational and rotational symmetry [311 113] , 
in attractive condensates [HIBS] . in low dimensions [^[50], for long-range interactions [51], 
and in metastable situations [521 [53]. In optical lattice systems, unusual depletion [51] and 
excitation frequencies measured in the presence of an external harmonic potential [HS] could 
not be explained within Bogoliubov theory. In all those situations, where the ground state is 
fragmented, the standard BdG approach for calculating excitation spectra is not applicable. 

A theory capable of describing statics of fragmentation phenomena is the best-mean field 
(BMP) method [37J. This method is based on a variational framework and a general mean- 
field ansatz for the state. It has been successfully apphed to evaluate the ground-state 
fragmentation of BECs in double-trap potentials [IQ] and allowed to identify stable frag- 
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merited excited states in repulsive condensates with energies below the GP self-consistent 
excited states [33] ■ The pathway from condensation to fermionization, passing a variety of 
fragmented states, has been demonstrated in [3S]. Similarly a variety of new Mott-insulator 
phases has been found for optical lattices |56j. This method has been extended to bosonic 
mixtures, showing interesting demixing scenarios [57]. A dynamical theory based on a simi- 
lar ansatz as the BMF is provided by the time- dependent multi- orbital mean field (TDMF) 
[58] method. For example, it allowed to predict interaction-induced self-interference fringes 



In this paper, we will generalize the standard BdG equations originally derived for con- 
densed (or simple [13]) BECs to the case of fragmented BECs. In particular, linear response 
of the TDMF will provide a tool for studying excitation spectra and properties of the ex- 
citations of fragmented condensates. The derived equations are general and can thus be 
employed for the calculation of spectra of systems with an arbitrary degree of fragmenta- 
tion. 

As an application of our response theory for fragmented BECs we will investigate Bose- 
Einstein condensates in symmetric and asymmetric double- well potentials. Those are pro- 
totypical systems which are best suited to test the theory, because they show fragmentation 
even for moderate interactions. Typically, double-well potentials have been studied either 
with a classical two-mode description [601 El] ('Josephson physics'), or within a two-mode 
Bose- Hubbard model [HI |62]. Other studies suggested that more than one state per well 
has to be taken into account [35l [63]. Recent dynamical studies [6H468] involve the mul- 
ticonfigurational time- dependent Hartree for Bosons (MCTDHB) method [361 EH]. A few 
excited states have been calculated in [301 [TO] using self-consistent methods. Bogoliubov 
approximation has been applied to few-site or lattice models [TH [72], and to double- well 
BECs [73]. Linear response studies beyond BdG have been performed so far with the sine- 
Gordon model [71] or within Gutzwiller-approximation [751477] . However, a thorough study 
of excitations is still missing and will be provided in this work. In particular we compare 
excitation energies and the density response of condensed and two-fold fragmented states in 
double-well potentials. 

The structure of the paper is as follows. First, in Sec. |TT] we introduce the underlying 
theoretical tools and models such as the standard BdG equations in Sec. II A and the TDMF 



method in Sec. II B In Sec. HI we present the derivation of the linear response theory for 



fragmented BECs, discuss its properties and provide expressions for important observables. 
Thereafter, in Sec. |IV]we apply the derived equations to calculate excitation spectra of BECs 
in double- well potentials. First, we discuss the possible structures of the ground states in 
this system, and then present our linear response results for symmetric and asymmetric 
double- well potentials in Sees. IV A and IV B respectively. Finally we summarize in Sec. [V] 
our findings and draw conclusions. There are also three appendices. Appendix |A] deals with 
some algebraic subtleties of the linear response equations. Special cases of linear response 
relevant for the application part are discussed in Appendix [B| where we also give explicitly 
the corresponding response matrices. In Appendix [C] we compare qualitatively the linear 
response of a two-fold fragmented system to that of a two-species system of distinct BECs. 
Finally, in Appendix [D] we derive the linear response of the two-site Bose-Hubbard model, 
which allows to quantify the importance of hopping excitations. 



II. THEORETICAL CONCEPTS 

A system of N interacting atoms in an external potential is described by the many-body 
Hamiltonian [T21IT3] : 

TV N 

^ = 5^Mr,) + Ao ^(r.-r,), (1) 

i=l i>j=l 

with the single-particle Hamiltoniai|^ 

h{r) = - VV2 + V{r) . (2) 

is the total number of atoms. The first term on the right-hand side of Eq. ([T]) describes 
kinetic and potential energy. The second term accounts for atom-atom interactions with 
interaction parameter Aq, which is proportional to the s-wave scattering length [121 113] . 
We use the commonly employed delta potential, but stress that the following formulas and 
derivations do not rely on the type of interaction potential. 

We work in dimensionless units where the energv is measured in terms of . We choose h = 1, set the 
mass of a *^Rb atom to one, and the unit of length to L ^ l/im. This gives units of energy and time 
E = h ■ 116.26 Hz and T = 1.37 ms, respectively. 
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A. Linear response of the Gross-Pitaevskii equation 



The standard method for calculating excitation spectra of interacting bosons at zero 
temperature is solving the Bogoliubov-de Gennes equations They can be derived 

as the linear response of the GP equation to an external time oscillating potential [I5l [78] . 
The linearized equations of motion are equivalent to the equations obtained when treating 
the many-body Hamiltonian, Eq. ([T]), in Bogoliubov approximation [iTl [79], [80], or in the 
random-phase approximation (RPA) [81j. All those approaches assume that the system is 
essentially condensed, i.e., one eigenvalue dominates the one-body reduced density matrix 
of the system [S2]. The non-condensate fraction has to be much smaller than unity in the 
Bogoliubov or RPA treatments. 

We shortly sketch here the derivation of the BdG equations [78]. The GP equation, which 
assumes that all atoms reside in a single orbital, reads 

# = ^GP0, HGP = h + \\(t>W (3) 

with the interaction strength A = Ao(A^ — 1). A small time-dependent periodic perturbation 
of the external potential, h{Y) — )■ h{Y) + 6h{r,t), can be written generally as: 

5Mr,t) = /+(r)e--* + r(r)e^'^S (4) 

with the probe frequency oo and the amplitudes realj^ Making the ansatz 

VN4>{r, t) = e"^^* \VN^\r) + u{r)e~'''' + t;*(r)eH , (5) 

as an expansion around the solution of the static GP equation (f)^{r) (with chemical potential 
fi) and with small amplitudes \u) and \v), one arrives at the equation 







1 1"> ) 





lCs^-^)\-\-\ f''^'"^\. (6) 



'BdG 



(7) 



The linear response matrix reads 

^ We note that without a perturbation, the fohowing procedure amounts to hnearizing the equations of 
motion Eq. (|3|. However, in order to show that the spectrum defined by the hnearized equations indeed 
corresponds to the frequencies of the excitation energies, we derive expUcitly the response to a small 
perturbation. The exact shapes of the perturbations /'''(r) and /~(r) do not influence the linear response 
spectrum. Rather the pole strength of the various excitations in the perturbed orbitals is affected. 
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Eq. ^ with the right hand side equal to zero is referred to as Bogohubov-de Gennes 
equations. They determine the response frequencies Uk, and also the response amplitudes 
d'u'^), \v''))'^. Using them we can solve the linear response Eq. ^ for {\u), with respect 
to a given perturbation. Inserted into the ansatz, Eq. ([s]), we obtain finally 

0(r,t) = e-'>^' |0°(r) + ^ E hu\r)e-'^' + 7>''*(r)e^"*] /{co - , (8) 

with the response weights (or pole strengths)]^ 

-fk = VN J rfr[M'='*(r)/+(r)0°(r) +t;'='*(r)r(r)0°'*(r)] . (9) 

In deriving the linear response theory for fragmented BECs we will keep the above nomen- 
clature as much as possible. 

B. Time-dependent multi-orbital mean field 

The Hilbert space of the many-body Schrodinger equation with Hamiltonian given in 
Eq. ([T]) is huge for the atom numbers one is interested in and which are typically used 
in experiments (say > 100). The following variational framework provides an efficient 
method to numerically solve this equation [371 EH]- The starting point is an ansatz for an 
arbitrarily fragmented mean-field state in terms of time-dependent orbitals: 

*(i"i, Yn, t) = 50i(ri, t)...0i(r„^, )f:)02( t)...(j)M{rN,t) . (10) 

Here, we put ni atoms into orbital 1, n2 atoms into orbital 2, and um into orbital M, with 
"^fLi ni = N . S is the symmetrization operator for bosons. Obviously, a GP state, where 



all atoms occupy one and the same orbital, is a special case of Eq. (10), i.e., for M = 1. We 
are thus about to generalize the GP equation. 



The ansatz Eq. ( 10 ) is now used to formulate an action functional 

M 



S = j rft|(vi/|^-,^|x[/)_^n,/i,,(t)[(0,|0,)-5,,]| , (11) 



The response diverges at the resonance frequencies, which seems to be unphysicaL When performing 
a time-dependent simulation, however, the response is damped due to effects beyond the hnear regime, 



which reestabhshes the physicaUy expected behavior |78) . 
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where fiij (t) are Lagrange multipliers which ensure the orthonormality of the time-dependent 



orbitals 0i(r,t). The first expression under the integral in Eq. (11) is evaluated to be 

rii hii - I i— I + \n—^ Wan + > XnuM ijij 



d 



i=l 



dt 



- + -^0^^ — Wan + 2J K^jWi. 



(12) 



with the matrix elements 

hii = j (l)*{r,t)h{r,t)(l)i{r,t)dT . 

and 



.d_ 
dt 



i / 0*(r,t)</>i(r,t)(ir , 



Wan = I \4>i{r,t)\^\^,{r,t)\^dr. 



(13) 



(14) 



The variational principle now requires the stationarity of the action, Eq. (11), with respect 
to the orbitals: 

0, i = l,...,M. (15) 



From this we obtain after some algebra the TDMF equations for a chosen occupation 

M 

Pi\4>i) = P /i + Ao(n,-l)|0,|2 + ^2Aon,|0,-| 



i), ^ = 1,...,M. 



(16) 



The projector P = 1 — X^fli \4>s){4>s\-, resulting from the Lagrange multipliers, keeps the 
orbitals orthonormal throughout the time propagation. The energy of the time-dependent 



mean-field states as obtained from Eq. ( 16 ) is conserved whenever the Lagrange multiplier; 



are hermitian, i.e., riijiijit) = njfi*^{t), or alternatively if {(pi\4>j) = 0, = 1, ...,M). Since 
there is no rigorous proof for the hermicity of fiij (t) throughout the propagation in time, one 
enforces the latter condition as an additional constraint. As a consequence, the projector on 



the left-hand side of Eq. (16) can be omitted. This has also the positive effect to simplify 



those integro-differential equations. With this we arrive at the final form of the TDMF 
equations: 

M 



t\(f),) = P 



h + Xo{n, - l)|0ip + ^2Aon,| 



^), ^ = 1,...,M. 



(17) 



^ Note that we use a different definition of tlie Lagrange multipliers as compared to Ref. [58] . The convention 
which is used there can be obtained by replacing /iij(t) — > iiij{t)/ni. 
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In this work we find that the hnear response of both the full form Eq. (16) and the working 



equations [see Eq. (17)] give rise to the same excitation energies and response amplitudes, 
and, therefore, to the same perturbed orbitals (f)i{r,t). Algebraic subtleties of the linear 
response of the full form are discussed in Appendix |Aj 

III. LINEAR RESPONSE THEORY FOR FRAGMENTED BOSE-EINSTEIN CON- 
DENSATES 

A. Derivation 

We now derive the linear response theory for fragmented BECs. For this purpose we add 



a time-dependent perturbation to the TDMF equations [see Eq. (17)], in a way as it has 
been done in Eq. ^ for the GP equation. The corresponding equations of motion can be 
written in a compact form as {i = 1, M): 

(z^-^g-^ 10.) -J2^'^M<pJ) = -mm , (18) 

where 

M 

Zi = h + Xoiui - l)|0i|2 + ^2Aon,|0,f (19) 

and 

/ii,(t) = (0,|Z, + 5/i(t)|0,). (20) 
The Lagrange multipliers ^ij(t) account for the orthonormality of the orbitals (pi, analogously 



as the projector P in the TDMF equations [see Eq. (17)]. Keeping fiij{t) explicitly helps 
to identify throughout the derivation the terms which originate from the orthogonalization 
with respect to the orbitals. We expand the perturbed orbitals around stationary solutions 
as 

0i(r,t)^0O(r) + <50,(r,t). (21) 
By keeping terms up to first order in 6(f)i{r,t) and /^(r) we obtain: 



/ f) X M 

f Z° -Z- + Ao(n, - m + ^2Aon,0f0°|<50,) + Ao(n. - iMfm) 

M M 

+ ^2Aon,-0>°|50*) - J2 Kl^*^.-) + ^f^^M<P')] = ■ (22) 



The zeroth-order contribution Zf is defined as in Eq. (19) but with 



\ Without the 



perturbation, i.e., 6h = 0, Eq. (18) is solved by the time-independent orbitals 



which 



are solutions of the best-mean field equations [37|. Those equations describe the stationary- 
states of the TDMFB 

M 

" (23) 



In many cases linear response is performed for ground-state orbitals, but 0^(r) could be 
excited stationary orbitals as well. The Lagrange multipliers to zeroth order are given as 



Hij = {(j)^\Zf\(j)^). The perturbed Lagrange multipliers are evaluated to be 



6nij{t) = 6 



Z, + 6h{t)) 10,) 



(24) 



M 



1=1 

In order to arrive at the first term we used that the unperturbed orbitals ^^(r) fulfill the 
best-mean field equations [see Eq. (23)]. Essentially, the matrix elements 6fiij{t) lead to the 



same projectors P as in TDMF, acting on Eq. (22). This is directly obvious for all terms 
of 6fiij{t) except for the first one. Using integration by parts and exchanging the indices of 
the summations, we can rewrite the first term of the sum YlfLi 



M 

j]/i°(50,i0?)i0;) 



M 



(25) 



We see that it corresponds to a projector on the term of Eq. (22) which is proportional to 
the Lagrange multipliers With this we find 



P 



Z^ + Xoin,-lM 



0|2 



M M 



M 



+Ao(n, - iMrm) + J]2Aon,0°0°|50*) 



-PSh{tM). (26) 



This equation has, similarly to the TDMF equations [see Eq. (17)], projectors on all terms 
except for the time derivative. 



^ We note that due to the presence of the Lagrange muhiphers, the stationary solutions of TDMF carry no 

time-dependent phase factors as they do in Eq. ([5|. 
^ Strictly speaking, the best-mean field is defined as the optimal orbitals at the energetically optimal 



occupation. 
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By using the ansatz 

^Mi{r, t) = ^/.(r)e-*'^* + <(r)e^"* (27) 



for the time-dependent perturbation to the orbitals {u is the probe frequency), and by 
equating hke powers of e^*"^*, we obtain the hnear response system 

where are the real amphtudes of the external perturbation. We switched to a vector no- 
tation in order to have a compact representation of the multi-orbital problem. For example, 
we denote the vector of stationary orbitals, multiplied by the square root of the population 
rii for each orbital, as = {\y/n^(f)^), l^yn^cp^), Ia/'t-a/^m))"^- ^ the linear response 



(29) 



matrix, with 

' Z°-H° + A B 



-B* -(Z°-^0'*) - A* 

Here, is a diagonal matrix containing the Z^. We group the matrix elements of the 
diagonal contributions originating from atom-atom interactions in A, as well as the off- 
diagonal ones in B. They are given as 



Ao(n,-l)|0Op, t = j 



i.,J^»<--^''^°'^' (30) 

[ 2Ao ^ 7^ J 

The projector matrix contains twice as many projectors as the number of orbitals M = 
1,...,2M) 

P , for i = j<M 

P* , for i=j>M , (31) 
0, {^^J) 



v., 



where P* = 1 — Yl!lLi l<^s)(</'sl- -^Y acting with (1 — P) on the linear response system (28), 
we find that for |ci;| > the solution is orthogonal to the stationary orbitals 0°(r), i.e.. 



(32) 







<\n)\ 
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This allows us to add an additional projector in Eq. (28) (i.e., replacing 'PC — )■ T'C'P). In 
order to find the excitation energies in Eq. (28) one has to solve the eigenvalue problem 





VCV ' = ' ' \ . (33) 



Most importantly and as a consequence of the projectors, for \uj^\ > each component of 
the eigenvectors (u'^,v'^)^ is orthogonal to the stationary orbitals 0°. 



We call the linear response Eq. (33), together with the linear response matrix given in 
Eq. (29), LR-BMF. The special case M = 1 is referred to as LR-GP, which is the linear 
response of the number-conserving version of the GP equation (for the differences to the 
linear response matrix of BdG see Appendix [B|. 



In order to find the orthonormalization relations of the response amplitudes for non-zero 



eigenvalues of Eq. (33), we study the symmetries of T'C/P similar as in Ref. [HQ]- First, a 



time-reversal spin-fiip-like symmetry: 

SiP£PSi = , (34) 

where the matrix = (5jj_A/ + 5i_A/j («, j = 1, 2M) permutes the i-th and the M+i-th 
raws, just as the first Pauli matrix cti = ( ? J) does for M = 1. Further, we have 

SsPrPSg = {VCV)^ , (35) 



where the matrix 




r , , ~. , , for i, ?' < M , , 

P3k = <i ' . (36) 

, for i,j>M 



For the case M = 1 this reduces to the third Pauli matrix = (q ^i)- From Eq. (34) 
we learn that, whenever (ju'''), jv'^))-^ is an eigenvector of 'PjC.'P with eigenvalue u'', then 
dv'^'*), lu'^'*))"^ is an eigenvector with eigenvalue —{u^)*. From Eq. (35) we find that 



(|u'^), — jv'^))"^ is an eigenvector of {"PCP)^ with eigenvalue to'', which allows us to construct 
the adjoint basis. From those considerations follow the biorthonormalization relations 

(vV*^''*)-(uV''*) = 0. (37) 
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It is obvious that all eigenvectors of 'PCV for = lie in the space spanned by the 
stationary orbitals of the unperturbed problem. The completeness relation then reads 



1 




M 



(Hi.o) + E 



({v'-''|,-{u'-'|) . 




(O.Kl)+E 



fc>0 




((u1,-(v^|) 



(38) 



(lu*^), jv'^))"^ are the eigenvectors with strictly positive eigenvalues uj^ > oj^ The i-th element 
of |u°) (|v°)) is equal to 0° (0°'*) = 1,...,M). All other elements of |u°) and |v°) 
vanish. 



Now we solve Eq. (28) by expanding the response vectors as well as the perturbation 
with the eigenvectors of CP orthogonal to the stationary orbitals 0° (r). The ansatz for 
the response amphtudes then reads 





(39) 



and for the perturbation 



-r(r)|0°) 




(40) 



Now Cfc and 7^ have to be determined. Substituting Eqs. (39) and (40) into Eq. (28), we 
obtain 

(41) 

where is defined in Eq. (|33|). From comparing coefficients in Eq. (|41|) we get an expression 





for the coefficients Cfc. Inserted in Eq. (39 ) this leads to a solution for the response amplitudes 
of the form 




E 



Ik 



U — UJk 




(42) 



Reinserting the amplitudes into the ansatz for the orbitals, Eqs. (|21|) and (|27|), we arrive at 

(43) 



the final solution for the time-dependent orbitals in linear response {i = 1, M): 

1 

/n. 



0,(r,t) = 0O(r) + 5^ 



/{UJ - LOk) . 



^ We found numerically real response frequencies u'' when starting from stationary orbitals. 



13 



Thus the orbitals, and with them the density, show the largest response at frequencies Uk- 
Moreover, the response for a fixed frequency is not necessarily equally strong for all the 
orbitals. This is because the components of the response amplitudes and Vj are not 



normalized, but rather the whole amplitude vector [see Eq. (37)]. The response weights, 
which quantify the intensity of the response, are given as 

M 

.fci*('„N -f+('„N J,0/'„N I _.A:,*/„\ ('„N J.0,*/ 



Ik 



E 



dr 



«f (r)/+(r)0;(r)+t;f (r)/-(r) 



(44) 



Similarly as the orbitals [Eq. (43)], it is dominated by the largest components of the response 
amplitudes. 



B. Density oscillations 

When probing the linear response through a time-dependent perturbation, an observable 
quantity is the oscillation of the density related to a given excitation [25j. From the orbitals' 



response, Eq. (43), we can calculate the time-dependent density for a probe frequency u and 
probing fields 

M M 

p(r, t) = Y, n^lMr, t) P ^ E ^^\^'(^)\' + 2 E Ap'(r) cos (cot) . (45) 

The density shows the largest response at the linear response resonance frequencies. For 
simplicity, we neglect here the typically very small imaginary parts of the response ampli- 
tudes, and assume real stationary orbitals 0° (r). We then obtain for the oscillatory part of 
the real space density 

M 
i=l 

The density in momentum space provides information about coherence properties of the 
system. For example, when two initially spatially separated parts of a BEC interfere in 
time-of-flight experiments, the density in momentum space describes approximately the 
interference pattern which is obtained on average. For a coherent BEC the fringe contrast 
is high, whereas it is zero for a two-fold fragmented BEC [83]. We note that in general 
interactions have to be taken into account during expansion and interference of the BECs, 
leading to interference fringes even for two independent BECs [59]. The Fourier transformed 
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orbitals and amplitudes, which we denote by </'°(p), u^^p) and (p), respectively, are in 
general complex. Therefore, the density oscillates at fixed amplitude but with a momentum 
dependent phase shift a'^(p): 

M M 

p{p,t) = J2nmP,t)\' ^ 5^n,|0°(p)r + 2 J]^^|Ap^(p)|cos [oot-a'ip)] , 
a^p) = arctan{lmAp'=(p)/ReAp'=(p)} . (47) 
The momentum-space density oscillations are given by 

M 

A-p'^ip) = J2 V^0-'*(P) K(P) + (P)] • (48) 



i=l 



Note that in Eq. (47) the modulus of the momentum-space density oscillations, |Ap (p)|, 
appears. This modulus can be measured in experiments as the maximal value of the density 
at each momentum p. 

The density response in position space can for some special cases be directly connected to 
the response weights. For a real periodic driving, which can be translated to /(r) = f^{r) = 



f (r), we can write Eq. (44) alternatively as 

7fc = /" drf{r)Ap\r) . (49) 



For this case the response of any observable is proportional to the density oscillations. 

IV. APPLICATION TO BOSE-EINSTEIN CONDENSATES IN DOUBLE- WELL 
POTENTIALS 

Before presenting our linear response studies of BECs in one-dimensional symmetric and 
asymmetric double-well potentials, we will briefiy discuss the structure of the (possibly 
fragmented) mean-field states in double-well potentials which are lowest in energy. 



Within a mean- field treatment, related to the ansatz of Eq. ( 10 ), the ground state in such 
a trap is either condensed or two-fold fragmented, depending on the barrier height and the 
interaction strength [HI]. The many-body wave function for the condensed state reads 



^(a;i,...,a;;v) = n,^i0°(x,), (50) 
whereas for the two- fold fragmented state with degree of fragmentation n/N it is given by 

v^(a;i,...,a;^) =5nr=i0°(a;Onf=,+,0°,(x,). (51) 
15 



Here, orbital (f>1{x) [4>%{x)] is localized in the left [right] well. The energies of those states 
are given by 



(52) 



and 



E 



M=2 



i=L,R 



+ 2\onLnR / \(l)l{x)\^\(l)l{x)\''dx^ 



(53) 



respectively. Above a critical barrier height, a fragmented state [Eq. (51) with n 7^ 0] 



becomes favorable in energy over a condensed one, Eq. (50). The same thing happens 



when the inter-particle interaction strength Aq exceeds a critical value. Typically, these 
transition points shift with atom number N (at fixed AqA^) to higher barrier and/or stronger 
interaction strengths. Importantly, even when the condensed state is lower in energy than 
the fragmented one, above a critical interaction strength the latter can be considered a stable 
excited state [53], [H]. It is typically slightly higher in energy than the condensed state, and 
is separated from it by an energy barrier. For example, in a symmetric double-well both the 
condensed and 50-50 left-right fragmented states are local minima with respect to a change 
in the critical occupation. In Fig. [T] we depict schematically these states and their energies 
for symmetric double-wells. 

In principle, besides the orbitals' excitations described by LR-BMF, there are excitations 
consisting of the redistribution of atoms between the orbitals ('hopping excitations'). Such 
processes can most easily be described within a two-site Bose-Hubbard (BH) model [4lJ . In 
Appendix D we derive the linear response of Bose-Hubbard ('LR-BH'), i.e., the response to a 
time- dependent potential as in Eq. (|4]). The excitation energies coincide as can be expected 
with the eigenenergies of the BH model. We study the LR-BH response weights and see that 
the dominant physical processes in a double-well potential with a time-oscillating potential 
perturbation are the orbitals' (or spatial) excitations. 

We note that the BMF and TDMF are mean field methods, and thus generally offer 
qualitative descriptions of BECs in double- well potentials [381 EEl [59] . In order to capture 
effects beyond mean field, one has to employ a description where the ground state of BECs 
in a double-well potential is neither completely condensed nor two-fold fragmented. In this 
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FIG. 1. (Color online) Schematic comparison between condensed, fragmented and stable excited 
fragmented states in a double- well potential. The single orbital of a simple BEC is delocalized 

(left chart). For a fragmented BEC the orbitals corresponding to different fragments [(^^{x) and 
0jj(x)] are localized (right charts). For a given barrier height, the stable fragmented state (lower 
right chart) is lower in energy than the condensed one up to some (possibly very high) atom 
number. Then, the condensed state becomes the lowest in energy. However, an energetically close 
stable excited fragmented state (upper right chart) typically exists. It is separated by an energy 
barrier from the condensed one. 

case the off-diagonal elements of the one-body reduced density operator (in the left-right 
basis) starts to play a crucial role The multiconfigurational Hartree for Bosons (MCHB) 
method [30J and its time- dependent variant, the multiconfigurational time- dependent Hartree 
for Bosons (MCTDHB) method [361 [69], offer full many-body descriptions. However, those 
methods are numerically much more demanding and are thus generally restricted to systems 
with smaller atom numbers and/or weaker interactions. 

In the following we study ultra-cold bosons in a one-dimensional double-well potential 
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parametrized as follows: 



V{x) = b/2 ■ cos {^x) + + a ■ X , (54) 



with the barrier height 6, harmonic oscillator frequency Uho determining the overall harmonic 
confinement, and asymmetry a. We present in the following excitation spectra of fragmented 



states, which originate from the derived LR-BMF response matrix Eq. (29) for the special 
case M = 2 (the linear response matrix is given in Appendix |B]) . To this we compare the 
response of condensed states, obtained from the number-conserving version of the BdG 
equations, i.e., LR-GP (those equations are discussed in detail in Appendix [B]) . Throughout 
this work we choose an harmonic confinement with uoho = v^- We discuss different values 
of the barrier height 6, as well as interaction strengths AqA^. 



The TDMF equations [see Eq. (17)] are for large enough atom numbers (say, > 20) 
practically independent of the total atom number A^, as long as AqA^ is kept fixed. They 
rather depend solely on the relative occupations rii/N. Similar statements hold also for the 
linear response. In this work we use A^ = 100 throughout, but we stress that the excitation 
spectra and corresponding observables are almost the same for larger atom numbers. 

Our linear response studies are based on a small perturbation of ground states. The 
ground state orbital for a simple BEG is obtained from the stationary GP equation 

Hgp<P\x) = fi<l)\x) . (55) 

For two-fold fragmented BEGs the ground-state orbitals are calculated as the lowest in 
energy solution of the M = 2 best-mean field equations |37] : 



hix) + XoiriL - l)|0^(x)p + 2Aon^|0°«(x)p| 0° (x) = (x) + , 

Hx) + XoiuR - iMlix)]' + 2Aonz.|0° (x)rt 4(a;) = /iV?j(^) + /^L0l(^) • (56) 



Those ground-state orbitals 0°(a;), 0° (x) and (p%{x) are real. 



Technically, we determine the ground states by imaginary time propagation of Eqs. (55) 
and (56) until the energy has converged to the desired accuracy of 10~^^. We use per orbital 



a grid size of A^ = 251 points and a box of size 12. The kinetic part is solved utilizing 
Fast Fourier transform. Using more grid points and/or a larger box does not lead to any 
visible differences in our plots. Then, to determine the frequencies Uk we diagonalize the 



non-Hermitian linear response matrix Eq. (29) on the same grid (per response amplitude). 
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We concentrate on the positive branch of eigenvalues (we recall that there is a negative 
partner to each eigenvalue, see Section III A). 



A. Symmetric double-well 



We start with a symmetric double-well potential, which is given by Eq. (54) with zero 
asymmetry a = 0. We first discuss the dependence of the excitations on the interaction 
strength and choose a high barrier in Sec. |IV A 1 Even for weak interaction strengths. 



the response of LR-GP and LR-BMF in momentum space is qualitatively different. For 
larger interactions two types of excitations emerge within LR-BMF and become energetically 
separate. Thereafter we proceed to study how the response changes with barrier height in 
Sec. lIVA2[ 



1. High harrier 

For high barrier heights the orbitals practically vanish around x = as can be seen in the 
inset of Fig. [2] As a consequence, the left (right) orbital of BMF has a shape similar to the 
left (right) half of the GP equation. Moreover, all terms in the LR-BMF response matrix, 
Eq. (29), which are proportional to |0/(;(a;)| ■ |0°(3;)|; are very small. Hence, the eigenvalue 



problem of Eq. (|33|) can be written as 

'0/ 



P 
P 



0\2„,fc 



uu 



'0/ 



0,* \ k 0,* k , ~ / jO,* 



y 



u, 



L ; 



-UVr 



(57) 



Similar equations hold for u'^ and (with indices L and R interchanged), see Appendix [B 
for the full matrix. We defined here n = XqN/2, and approximated ~ A^ — 1. Most 
importantly, u'l and are governed by the same operator as and v^. 



h + 2n{\<l>1\' + m . (58) 
Note the difference of Z^'' to the TDMF-operators Z^, which carry the index L or R. Hence, 



the first term in each line of Eq. (57) describes a particle in the effective potential of two 



condensates, both carrying a factor of two due to exchange interactions [58]. The term 
proportional to the off-diagonal Lagrange multipliers and /x^jl is a coupling term to the 
other amplitude m^. The last terms on the left hand sides couple to f^. 
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FIG. 2. (Color online) Excitation spectra of a BEC in a symmetric double-well potential (a = 0) 
versus interaction strength with high barrier height (b = 20). The total atom number is = 100, 
and the overall harmonic confinement given by ujho = V2. We compare the linear response of 
LR-GP, shown by the blue solid lines, and LR-BMF, shown by the red dashed lines for direct, 
and orange dotted lines for swapped excitations. The excitations are grouped into pairs of lines 
with gerade-ungerade symmetry and marked with numbers. The swapped excitations of LR-BMF 
are marked with primed numbers. Inset: We plot the potential by the black solid line. The 
corresponding ground state orbital of GP for XqN = 10 is shown by the blue solid, and the left 
orbital of BMF by the red dashed line. All quantities are dimensionless. 

For weak interaction strengths, the effects of coupling of to and v'l (and similarly 
for L R) are negligible. Since is symmetric in L and R, it originates to delocalized 
response amplitudes which have either gerade or ungerade symmetry]^ LR-GP reduces to 
an equation similar to that for or of LR-BMF. Hence, the energies [see Fig. [2] and 
® We note that for smaller atom numbers, on the order of iV = 100, the difference between N and iV — 1 

leads to localized orbitals. However, in this regime the two basis sets (i.e., left-right or gerade-ungerade) 

lead to the same physics. 
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response amplitudes coincide. However, as we will show later, the momentum space density 
responses of LR-GP and LR-BMF differ strongly due to the different structures of the ground 



states, i.e., coherent or fragmented [see Eqs. (50) and (51), respectively]. 

For larger interaction strengths, we find that for LR-BMF two types of excitations become 
energetically separated. In particular, an excitation can be either to an orbital, which 
dominates in the same well ('direct' excitation), or in the other well ('swapped' excitation). 
We sketch the notion of swapped excitations in Fig. [3j While the orbitals (pL and (pR are 
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FIG. 3. (Color online) Sketch of the two types of excitations occurring for fragmented BECs in a 
double-well potential. Direct (swapped) excitations can be interpreted as the promotion of a boson 
from one orbital to an excited state which dominates in the same (other) well, e.g., from (pi to 
{u\). The transfer of atoms to direct (red solid line) costs less energy than to swapped (orange 
dashed lines) excitations. This is because a depletion, related to (green dashed-dotted line), 
reduces the energy of direct excitations. For the swapped excitations, v]^ « 0. Response weights 
of direct and swapped excitations of LR-BMF are proportional to the overlap between an orbital 
and the corresponding response amplitudes. Hopping excitations of LR-BH are proportional to the 
much smaller overlap of the two ground-state orbitals. 



localized and have very small overlaps, the w-amplitudes of LR-BMF are partly delocalized. 
We can understand the energetical sphtting between direct and swapped excitations as 



follows. In Eq. (57), the term which accounts for a coupling to becomes important 
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for larger interactions. The w-amplitudes describe a depletion of the true ground state 
of a condensate in which a few atoms occupy excited states [HI]. Furthermore, we find 
the f-amplitudes, and, hence, the depletion, to be local. They are nonzero only for direct 
excitations, see Fig. |4} The depletion thus lowers the energy of direct excitations as compared 
to swapped excitations where the energy is determined solely by exchange interactions [see 
red dashed and orange dotted lines in Fig. |2| respectively] . We conclude that we found a 
class of excitations in a fragmented system which do not appear at all in a condensed one. 

Accordingly, for a fragmented state the amplitudes are localized, in contrast to LR-GP. 
The response amplitudes have either gerade (g) or ungerade (u) symmetry, which we label 
as k = Ig, 1m, where k is the index of the excitation. We show the gerade ones u^^{x) and 
v^^{x) in Fig. (a) by the blue solid and dotted lines, respectively. The linear combination 
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FIG. 4. (Color online) Example of response amplitudes u'^ and v*^ for the same double-well system 
as in Fig. [2] and for Aq^V = 10. (a) The amplitudes u and v of LR-GP, shown by the blue solid and 
short-dashed lines, respectively, are symmetric functions. For LR-BMF, instead, they are localized 
functions. We show u]^ by the red dashed line, and v]f by the red dotted line, (b) Amplitudes of 
the swapped excitations Uj^ and of LR-BMF, shown by the orange dashed and dotted lines, 
respectively. Amplitudes of LR-BMF are magnified by a factor of \/2 for easier comparison. See 
text for more details. All quantities are dimensionless. 

of left and right LR-BMF amplitudes is either gerade or ungerade, and we thus use the same 
labeling as for the condensed state. We show in Fig. |4] (a) the first direct excitation of LR- 
BMF. In particular we plot the LR-BMF amplitudes u]^{x) and v^^{x), which have for a; < 
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the same shape as for LR-GP (except around x = 0). We further note that the response 
amphtudes in both condensed and fragmented cases have a node which ensures orthogonahty 
with respect to the ground-state orbitals. In Fig. |4] (b) we show the first swapped excitation 
of LR-BMF, i.e., the amphtudes u\^{x) and f^^(x). The response amphtude f^^(x) almost 
vanishes. 



The terms of Eqs. (57) proportional to the Lagrange multiphers ^Ij^ and /i/j^, as well as 



the terms proportional to the orbitals as ■ which we neglected in Eq. (57), induce 
another small shifts of energies. 

The response weights of swapped excitations are given by the overlap integrals of orbitals, 
response amphtudes, and perturbations: 

71' = / dx [4'*r<Pl + n];'7V?, + VL'*r<Pl + vn'*r<pi) ■ (59) 

Although it comprises a transfer of atoms to the other well, it is in general dominant over 
excitations involving a hopping of atoms between the orbitals. Within linear response of the 



Bose- Hubbard model ('LR-BH', for details we refer to Appendix D) and under a periodic 
potential perturbation, hopping can occur either directly through tunneling, or mediated 
by a time-dependent potential difference between the wells. The first process is completely 
irrelevant here since the orbitals (t)\{x) and (f)%{x) have small overlaps. The second process 
has finite response weights only for the first few hopping excitations, with energies well below 
the first excitation of LR-BMF. 

We move on to the study of excitations by means of an observable. Let us first discuss the 
limit of very weak interaction strengths by means of approximate formulas for the orbitals 
and response amplitudes. In particular, we denote with tp"'{x) {n = 0, 1, 2, ...) the normalized 
ground and n-th excited harmonic oscillator eigenfunction, centered around x = 0. We model 
the GP orbital and first excited w-amplitudes of LR-GP as (the f-amplitudes are negligible 
for weak interaction strengths): 



(f>\x)=[^lj\xL)+^\xR)] /V2, 

u'9,u^^^ ^ [4j\xL)Ti^\xR)] /V2, (60) 



with xl '■= X + d {xr := X — d). d is half the distance between the minima of the left 
and the right wells. The minus in front of iP^{xr) is needed in order to construct a gerade 
function out of two displaced ungerade functions iP^{xl) and 'iP^{xr). The BMF orbitals 
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and LR-BMF response amplitudes can be considered to be completely localized. For direct 
excitations we have 



0° (x) = , = ij\xn) , 

u'l\x)=ij\xL)/V2, uf\x) = Ti^\xn)/V2. (61) 
The normalization of the u-amplitudes follows from the orthonormalization relations in 



Eq. (37). The density oscillations of the direct excitations are obtained by plugging Eqs. (60) 



and (61) into Eq. (46): 



^ ^ = [^P\xl) + ^\xn)] [^P\xl) T ^\xn)] /2 



N 



BMP 



N 



[^lj''ixL)^\xL) T ^'ixn)^p\xn)] /2 . (62) 



When assuming that the overlap between displaced functions vanishes, it results that 
[Ap^9,«]^^ = [Api3,«]^^^^. Similarly, for the response weights holds big.Jcp = [lig,u]BMP- 
Moreover, the density response of the swapped excitations of LR-BMF and their response 
weights vanish. Thus, for large barriers and weak interaction strengths, the density in 
position space responds in exactly the same fashion for both condensed and fragmented 
states. 

But what if we proceed to momentum space? In this case, the ground state densities 
are qualitatively different: whereas for GP the density shows up a modulation due to the 
coherence between the bosons in the left and right well, the density of a fragmented state is 
simply a Gaussian. Using a similar notation as above, we denote the Fourier transformed 
ground and excited harmonic oscillator states as ip"'{p) (n = 0, 1, 2, ...). We remind that the 
Fourier transform of a translated function amounts to the Fourier transform of the original 
function times an additional phase factor, i.e., FT [ipn{x — d)] = e~^^'^FT [ipn{x)]. We then 



find for the condensed system as Fourier transforms of Eqs. (60) 



0°(p) = V2 COS {pd)ip%p) , 
u^^ip) = y/^i sin {pd)tjj^{p) , 

u^"(p) = V2 cos {pd)ip\p) , (63) 
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and for the fragmented one as Fourier transforms of Eqs. (61) 



u'r{p) = Ti^\p)e-'''' 



(64) 



Plugging Eqs. (63) and (64) into Eq. (48), we obtain for LR-GP the following density 



oscillations in momentum space 

m''ip)]Gp\ 



N 



|[Api«(p)] 



GP\ 



N 



w{p) '^/'^(p) sin (2p(i) 



2V^°(p) ^^{j)) [cosM] 



(65) 



which show up modulations of the phase with frequencies 2pd and pd, respectively. For the 
first direct excitation of LR-BMF, marked as 1, we have 



[Api'^p)] 



BMP 







[Apl'"(p)]BMF 



(66) 



Thus, the momentum-space density response vanishes for the gerade direct excitation and 
the ungerade direct one has only one node. Thus, even for a high harrier and weak interaction 
strengths, the momentum-space density oscillations of condensed and fragmented BECs are 
qualitatively different. 

For completeness we quote here also the momentum-space density response for the first 
pair of swapped excitations, marked as 1', although the corresponding response weights 
vanish for weak interaction strengths. It can be obtained by switching the signs in the 



exponents of the second and third quantities of Eqs. (64). Plugged into Eq. (48) this results 
in 

\[^~P^^{p)]bmf\ 



N 



|[Api«(p)] 



BMF\ 



N 



w{p) w{p)cos{2pd) 



(67) 



We find that the gerade-type density response of LR-GP and the swapped gerade one of 
LR-BMF are very similar and have the same period. This is in contrast to the ungerade-type 
response, where for the swapped excitations of LR-BMF the period doubles. 
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We now proceed to the density response at stronger interaction strengths. In Fig. [5] we 



plot the density oscillations in position space as defined in Eq. (46) for different excitations. 
The broad gray solid and dark dashed lines in panels (a) and (b) show the ground state 
densities for GP and BMF, respectively. They perfectly coincide for high barriers. In the 
left panels we show results for AqA^ = 10, and in the right ones for AqA^ = 20. It is interesting 
to note that the f -amplitudes have typically opposite signs as the w-amplitudes (see Fig. 
They become important for large interaction strengths, and lead to a damping of the density 
oscillations (compare left and rights panels in Fig. |5]). For the lowest excitation, marked 
as 1 [see Fig. |5] (a) and (b)], we find that the density oscillations of the condensed (blue 
solid lines) and fragmented states (red dashed lines) have almost the same shapes. We plot 
only the gerade density oscillations, but there exist counterparts of ungerade type as well. 
Whenever gerade and ungerade excitations are energetically degenerate, the moduli of their 
real space density oscillations are on top of each otherj^ 

For the swapped excitations of LR-BMF, the response amplitudes are partly delocalized 
[see Figs. |3]and|4] and thus the position-space density response is non-negligible even for the 
lowest swapped excitations at A = 20, marked as 1' [see Fig. [s] (b)]. Similarly, the response 
weights are non- negligible. 

In momentum space, the GP ground state, plotted by the gray line in Fig. |6] (a), shows 
an interference pattern refiecting the coherence between left and right condensates. This 
is completely different for BMF, plotted by the dashed black line, which describes two 
independent condensates. The momentum-space density response of LR-GP and LR-BMF 
for the excitations marked as 1 and 1' [see Fig. [6] (a,b,c)] are at larger interactions still 



qualitatively described very well by Eqs. (65), ([67]) and (66) (note that the behavior around 
p = is not captured well by the simple equations for the orbitals and amplitudes, and the 
gerade solutions do not vanish at this point). 

In the next higher group of excitations, marked as 2 and 2', also the direct excitations 
of LR-BMF deviate from the LR-GP energies, as can be seen in Fig. |2j For LR-GP we find 
a splitting between gerade and ungerade excitations. The excitation energies of LR-BMF 
lie between them. The position space density response of the swapped excitations becomes 
more sizable as for the lowest excitation, see Figs. IS] (d,e,f). Thus, swapped excitations 



The response at a degenerate frequency is the sum of both the gerade and ungerade contributions, mul- 



tiplied by the corresponding response weights Eq. ( 44 1 . 
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become important for excitations with energies of the order of the barrier height. Whenever 
the swapped excitations give rise to an appreciable position-space response, this signifies the 
transition from below the barrier to above the barrier excitations for a fragmented system, 
inasmuch as the hfting of the g-u degeneracy does for the response of a condensed system. 

For the excitations marked as 3 and 3', the amphtudes become now more and more 
delocahzed. This is a consequence of the fact that the higher the energy of the excitation, the 
less important is the barrier. While becoming similar in size, the direct excitations dominate 
in their own well, and the swapped ones in the respective other well. As a consequence, 
within BMF we find a position-space density response weaker by a factor of a/2 as compared 
to LR-GP, see Fig. [s] (e,f). For weak interaction strengths this can be understood replacing 



in Eqs. (61 ) 



= u'rix) = [^PHxl) T ^Hxn)] /2 . (68) 

From this we obtain [Ap^f'"]^j^,^^ = [Ap^f'"]^p /a/2, as well as [yig^ulsMF = [lig,u]Gp/ V^- 
Quahtatively, these formulas describe also the density response for stronger interaction 
strengths. 

To gain further insight, we plot in Fig. [7] both gerade and ungerade position-space density 
oscillations for LR-GP with AqA^ = 20, as shown by the blue solid and dashed-dotted lines, 
respectively. The shapes of the gerade and the ungerade excitations are different, reflecting 
the lifted energetical degeneracy of the excitations. For the fragmented system we have two 
pairs of degenerate excitations: direct and swapped ones. We plot the gerade direct density 
response by the red dashed line, and the ungerade swapped one by the orange dashed- 
double dotted line. For a simpler comparison between LR-GP and LR-BMF we compensate 
for the factor of a/2. We conclude that well above the barrier the direct excitations of LR- 
BMF become similar to the gerade ones of LR-GP, although they have different energies 
[see Fig. [2], and the swapped excitations become similar to the ungerade ones. Similar 
statements hold in momentum space, see Fig. [6] (g,h,i). 

For LR-GP there is in principle also an excitation where an atom is transferred to the 
ungerade solution of the GP equation. However, for the case of a high barrier the gerade and 
ungerade solutions of GP are almost degenerate, even for large interactions. Therefore this 
excitation has very small energy for all values of A. Since it cannot be distinguished from 
zero in the plot, we do not show it in Fig. [2j For the linear response of a BEG prepared in 
the ungerade GP state, the parity of the position space density oscillations changes. For the 
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momentum space density oscillations, the periods of gerade and ungerade would be reversed 
as compared to Eqs. (65). 



2. Low barrier 

with it Now we study what happens if we lower the barrier. In this case the orbitals of 
course penetrate the barrier more than at high barriers, see the inset of Fig. [8] Here we 
observe an excitation of LR-GP which lies energetically below all other excitations (labeled 
0). This excitation corresponds to the ungerade solution of GP. When decreasing the barrier 
height, the energy of excitation increases. This signifies, that for low barriers the ground 
state is condensed. 

However, also fragmented states in a trap with low barrier could be of relevance in 
experiments. For example, when one cools down a thermal gas in double-well potential, 
it is not clear whether one manages to condense the atoms really into the ground state of 
the system, or the system resides in stable fragmented excited states [85]. Moreover, when 
two initially independent condensates are slowly merged, the question arises, whether the 
system evolves adiabatically and becomes condensed, or stays fragmented. When analyzing 
excitation spectra of a fragmented BEG using LR-BMF, the excitation is absent, because 
there is no phase relation between the left and right condensates. Thus, the existence of an 
excitation to an ungerade mode can be used as a signature for coherence in the system and 
to characterize the state obtained by cooling into a double-well potential. 

For the higher excitations of LR-GP, we clearly see how they become degenerate as the 
barrier is raised. For example, around 6 ^ 18 the lowest pair of excitations, marked as 1, 
becomes degenerate (blue solid lines). For a fragmented ground state the situation is differ- 
ent: At small barriers (say b < 15), LR-BMF has twice as many non-degenerate excitations 
as compared to LR-GP. This splitting is due to the terms in the response matrix which 
are proportional to the overlap of the ground-state orbitals 0° 0r (see Appendix [b|) . 
Already around 6 ^ 15 states become pairwise degenerate. In contrast to LR-GP, the two 
branches of BMF excitations stay energetically well separated. 

The position-space density oscillations at low barrier 6 = 14 are shown in Fig. |9| Ex- 
citation as shown in (a) by the blue solid line is qualitatively similar to the GP ground 
state, but it has ungerade symmetry and a flattened top. The excitations marked as 1 and 
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1' are shown in (c,d) and (b), respectively. Most importantly, while the gerade excitations 
marked as Ig are similar for LR-GP and LR-BMF, the ungerade ones marked as lu are 
quantitatively different [see (c) and (d), respectively]. Moreover, the swapped excitations 
of a fragmented state are approximately as large as the direct ones. Also for the density 



oscillations in momentum space. Fig. [T0| we find strong differences for condensed and frag- 
mented systems. Hence, the lower the barrier and thus the larger the overlap of the left and 
right condensates, the more striking are the differences of LR-GP and LR-BMF. 



B. Asymmetric double-well 

After having applied our response theory to BECs in a symmetric double well potential, 
we now turn to a slightly asymmetric double- well. We use an asymmetry of a = 0.1 and 



choose a relatively high barrier height b = 20 [see Eq. (54)]. The corresponding potential is 



shown in the inset of Fig. 1 1 



Condensed and two-fold fragmented mean-field states compete for being lower in energy 
also in an asymmetric double well [HI]. The condensed one extends over both wells and 
dominates in the lower well. The fragmented BEG consists of a larger fragment in the lower, 
and a smaller one in the upper well. We show GP (blue solid line) and BMF (red dashed lines) 
orbitals for typical parameters in the inset. For the chosen interaction strength AqA^ = 10 
we find that the lowest in energy mean-field state is two-fold fragmented, except for very 
large atom numbers. In the latter case the condensed state is slightly lower in energy and 
the fragmented state becomes a stable excited state [HI]. The optimal occupation difference 

— which characterizes a stable fragmented ground or excited state depends in principle 
on N. However, it takes on practically the same values already for > 100. The ground 



state energy versus occupation n^, is shown by the filled gray line in Fig. [11} showing a 
minimum at ul/N ^ 0.55. The ground state densities of GP and BMF, which are shown in 
Fig. [121 perfectly coincide. 



Let us first discuss how the excitation frequencies depend on the BMF occupation {n^, N- 



hl) as shown in Fig. 11 Remarkably, the first excitation (marked as 1) is independent of the 
occupations. We can attribute this to the fact that the response amplitudes of the fragmented 
system are purely local in this case (not shown). They are, therefore, determined by the 
local ground state density, which is mostly independent of the occupations for the range 
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of riL as shown in the figure. In contrast, the swapped excitations (marked as 1') depend 
hnearly on rti and cross each other approximately at the optimal occupation ni/N = 0.55. 
This is due to the fact that the energy needed to excite an atom to the first excited state of 
the other well depends on the chemical potential difference fiiL — f^RR [see also the matrix in 



Eq. (B5) in Appendix |Bj. This quantity has been identified as the energy needed to transfer 
a boson from one well to the other (at large enough N) [81] , 

= /ULL - ^J'RR ■ (69j 

We observe that at the optimal occupation the transfer of bosons is suppressed and, therefore, 
IJ-LL = fJ'RR- Hence, in addition to the exchange interactions also the potential difference 
contributes to the energetical splitting of direct and swapped excitations. The higher lying 
direct excitations cross each other, similar to the swapped ones. This is because as soon as 
the amplitudes of the direct excitations become delocalized, the energy depends on fiiL—f^RR- 
We next discuss the energies and the density response at the optimal BMP occupation 



rii/N = 0.55 (marked by the dotted vertical line in Fig. 11). For the lowest excitation 
(marked as 1) we observe that the response frequencies for both condensed and fragmented 
systems (direct excitations) coincide and are doubly degenerate. The corresponding density 



oscillations in position space for the two solutions are shown in Fig. 12 (a) for LR-GP 
and in (b) for LR-BMF (see the solid and dashed lines, respectively). Most importantly, the 
density response of a simple BEG is delocalized, while it is strictly localized for a fragmented 
BEG. Also the response amplitudes are localized in the sense that either {\u\} , \v'l))'^ or 
(|m|.), |f^))^ are finite, respectively. We label the two degenerate frequencies as k = la, lb. 
The total response at this frequency is then given by the sum of two contributions [see 



Eq. (45)], proportional to 

p{x) ~ ^laAp^'ix) + ^ibAp'\x) . (70) 

For weak interaction strengths we can model the imbalanced ground state orbital and 
M-amplitudes for the condensed system as 



0°(x)= V^^P^xl) + VN- Ul^^Xr) /VN 



{xl) ± - riLtp' {xr)] /y/N . (71) 
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The following approximations are particularly good for small — nj^. 
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For the strictly localized orbitals and excitation amplitudes of the fragmented system we 
have 

0° (x) = , = i^^XR) , 



u 



< = 0, < = ^i(x«). (72) 

From this we arrive at the delocalized density oscillations and response weights for LR-GP 
(assuming that the overlaps of displaced functions vanish): 



= / dxf+{x) [nL^''{xL)tp\xL) ±{N- nL)4j''{xR)4j\xjt)] /TiV. (73) 
For LR-BMF the same quantities are localized and given as 



[^p"'{x)]^,,^ = V^^P'{xL)^P\xL), 
7iT^ = V^ [ dxt{x)i,\xL)i^\xL) , 



„ BMF 
lib 



VN-ul / dxf-'{x)ij\xR)^\xR) . (74) 



The total density response of LR-GP is thus according to Eq. (70) proportional to 



PGp{x) ~ 2 [{nL)^L^%XL)^\xL) + {N- nLfiR^%XR)^\xR)] /N, (75) 

where we defined ^l,r = J dxf^{x)ilj^{xL,R)'>p^{xL,R) [f~{,x) does not appear because the 
w-amplitudes are marginal for weak interaction strengths]. For LR-BMF we obtain 

Pbmf{x) ~ nL^L'ip^{xL)ip^{xL) + {N - nL)^Ril)^{xR)il?-{xR) . (76) 

Hence, in general the total response of the lowest excitation is different for LR-GP and 
LR-BMF even for weak interaction strengths. The difference in the densities is proportional 
to the imbalance ul — ur. It vanishes for symmetric occupations rii = N/2. In this case, 
the above equations boil down to the results for a symmetric double-well potential as given 
in Sec. WJ\ 

For the swapped excitations the amplitudes are finite in that well where the corresponding 
ground state orbital vanishes [orange dashed and dashed-dotted hues in Fig. [12] (b,d,f)]. 
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Similar to the direct excitations, the position-space density response of the swapped ones is 
strictly localized. However, the left and right swapped excitations have different energies, 



and the position-space density response for the lowest two of them, marked as 1', see Fig. [12 
(b) (orange dashed and dashed-dotted lines), is very small. 



We next discuss the density response in momentum space, see Fig. 13 (a-c). For LR-GP 
(blue solid lines) it is similar as in the symmetric double-well [compare with the respective 
results plotted in Fig. [6] (a-c)] . For LR-BMF, however, the momentum-space density oscilla- 
tions are very different as compared to the symmetric case - they have very little structure 
and only one node. This can be explained by the local nature of the excitations. For weak 
interaction strengths the momentum-space density oscillations of the direct as well as the 



swapped excitations are, following from Eq. (72), just given by 

\[Ap^{p)]bmf\ = V^L^\p) \i'\p)\ , for A; = la, 1'6, 

\[/^~p\p)]bmf\ = VN-nL^p\p) \i'\p)\ , for k = lb,l'a. (77) 

Hence, the phase factors which have been found in case of a symmetric double-well 



[Eqs. (66)], and which are due to an interference of the left and the right response, are 



absent here. 



For the higher excitations we find an energetical splitting of the lines, see Fig. 11 For 
LR-GP the density oscillations with an imbalance on the left (blue solid lines in Fig. 12) 
or right (blue dotted lines) occur at different frequencies. Similarly, the left (red solid and 
orange dashed lines) and right (red dotted and orange dashed-dotted lines) localized density 
oscillations of LR-BMF become energetically separated. We note that the shapes of the 
density oscillations are strongly affected by the contribution of the off-diagonal Lagrange 



multipliers piji and in the linear response equation [see Eq. (57)]. Hence, in contrast to 
a condensed state, the position-space density response of a fragmented state is purely local 
with generally different frequencies and density oscillations for the left and right response. 



V. SUMMARY AND CONCLUSIONS 



The Bogoliubov-de Gennes equations are the standard linear response equations for 
bosons. They are applicable for condensates where all atoms occupy only a single orbital. 
We presented the linear response theory for fragmented condensates, where the atoms are 
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allowed to be distributed over several orbitals. We derived the linear response equations for 
a small periodic perturbation of a stationary state. Since our linear response theory is based 
on the BMF and TDMF methods, we call the derived equations LR-BMF. Those allow us 
to determine excitation energies and response amplitudes of fragmented condensates with 
an arbitrary degree of fragmentation. 

We analyzed the properties of LR-BMF. Most notably, in the derived equations the re- 
sponse of each fragment is orthogonal to all the ground-state orbitals. This has vast conse- 
quences on the shape and the energies of the excitations. The linear response matrix defines 
a biorthogonal basis, which consists of the vector of response amphtudes related to all the 
fragments. The response of the fragments is coupled through the Lagrange multipliers, and 
whenever they overlap in space. The Lagrange multipliers of the BMF ground-state orbitals 
enter the linear response matrix and account for the relative energies of the stationary or- 
bitals. We give expressions for the density oscillations in real and momentum space which 
arise due to a resonant perturbation. They are given as sums of the contributions due to 
the response of each fragment, weighted by the square root of the orbitals' occupations. 

As applications, we investigated Bose-Einstein condensates in symmetric as well as 
asymmetric double- well potentials. We compare results of the LR-GP (i.e., Bogoliubov-de 
Gennes) and LR-BMF theories. Our numerical results demonstrate that the responses of 
a fragmented and a condensed system show striking differences. In particular, fragmented 
BECs possess a class of swapped excitations which do not exist in condensed systems. 
They are characterized by a response which is dominantly in the respective neighboring 
well. These excitations signify the transition from below the barrier to above the barrier 
excitations. The density response in momentum space has been found to be qualitatively 
different, even in situations where the response energies of the two theories numerically 
coincide. At low barrier heights, an excitation to the ungerade state of GP exists within 
LR-GP, but it is absent within LR-BMF. Thus it can be used as a signature of coherence. 
For fragmented BECs in asymmetric double-well potentials we found a locahzed density 
response (i.e., finite in either one or the other well), as well as an energy splitting between 
left and right response. This is in stark contrast to the response of a condensate which is 
not fragmented. 

We conclude that, for a proper analysis of the response of even one-body observables like 
density, it is crucial to take into account the many-body structure of the underlying state. 
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In view of the rich physics which has been found using the standard Bogohubov-de Gennes 
equations, our generahzation of this very successful theory to fragmented BECs offers even 
more rich prospects for understanding excitations of cold atom systems in general. The vast 
differences between the response of condensed and fragmented systems will provide a way 
to distinguish condensed and fragmented states by linear response experiments. 
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Appendix A: Linear response of TDMF without constraint 



We discuss here shortly the derivation of linear response of the full form of the TDMF 



equations [see Eq. (16)]. The projector on the left hand side of Eq. (16) translates to a 



projector on the term proportional to u) in the linear response equations [see Eq. (28) 
As a consequence, the response amplitudes (|u), |v))"'" are not necessarily orthogonal to the 



ground-state orbitals 0° [i.e., Eq. (32) does not hold]. Thus, the question arises if in addition 



to the orthogonal eigenvectors as defined in Eq. (33), also the ground-state orbitals have 



to be included into the ansatz of the response amplitudes in Eq. (39). Since we expanded 



around stationary BMF (ground state) orbitals 0°, we assume that they are recovered if 
the perturbation is zero (i.e., = = 0). Thus, the response amplitudes (|u), |v))^ can 
contain solely those ground state orbital contributions, which lead to trivial time- dependent 
phases on the orbitals. We note that those are determined from BMF only modulo a phase. 



Thus, also when we start with the full form of the TDMF equations [see Eq. (16)], the 
frequencies Uk (excitation spectra), the response amplitudes |u'^) and jv'^), as well as the 
perturbed orbitals 4>k are the same as for the linear response of the TDMF working equations 
[see Eq. ^l^ . 
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Appendix B: Special cases of linear response: M = I and M = 2 



Linear response of a condensed state (M = 1) 

For M = 1 we recover the results for the excitation spectra as obtained from the number- 



conserving Bogoliubov theory of Ref. [80], with the linear response matrix from Eq. (29) 



'M=l 




(Bl) 



We stress that in our derivation the orthogonality of the response amplitudes {\u), \v))'^ to 
the ground-state orbitals (p^ is obtained naturally from the derivation, in contrast to Ref. [50] 
where it is an assumption. 

Interestingly, the standard and the number-conserving BdG equations can be considered 
as the linear response of different forms of the GP equation, which deviate by a global phase 
on 0(r, t) [58]. Such a phase has no physical meaning, and, therefore, the three forms of the 
GP equation are equivalent and predict the same physics. The standard form, with response 
matrix given by Eq. ([T]), is obtained from the usual form of the GP equation: 

z0 = Hgp^ . (B2) 



If we consider the TDMF equations [see Eq. (17)] for the case of M = 1, we obtain a GP 



equation which is similar to that obtained from a number-conserving approach (see Ref. [80]): 

i0 = PHGP(f) . (B3) 



The linear response of this equation leads exactly to the linear response matrix of Eq. (Bl). 
Both the standard and the number-conserving BdG equations lead to the same response 
frequencies u^. The response amplitudes in the number-conserving formalism are obtained 
from those of the standard form by projecting into the subspace orthogonal to the ground 
state orbital 0°(r) with the projector P [SD]. The difference in the perturbed orbitals 0(r, t) 
is then just a trivial phase. However, the zero eigenvectors of both forms differ from each 
other. In particular, in the standard, number-non-conserving linear response there is one 
missing zero eigenvector, which is supposed to lead to a divergence of quantum fluctuations 
in time [80l[86]. 
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The full form of the TDMF equations [see Eq. (16)], including the projector on the left 
hand side, reads 

Pt^ = pHGP<p, (B4) 



having the same linear response as Eq. (B3), apart from contributions of the ground state 



orbital to the solution of the perturbed orbital [Eq. (|8])], which corresponds to a physically 
irrelevant global phase. 



Linear response of a two-fold fragmented state (M = 2) 



We explicitly write down the linear response matrix which is used in the application part. 



Sec. IV, as a special case (M = 2) of the general formula, Eq. (29). Since the orbitals of a 
two-fold fragmented BEC in a double-well potential are localized, we use as orbitals' indices 
left 'L' and right 'R'. We obtain as the linear response matrix Cm=2 = ^m=2^- 



-'M=2 — 



12 



Tin ^ LR 



(B5) 



n J.0,* iO,* 

-2n0^ 0L 



' riR 0,* 



rV + f^RR ) 



where we use the notations rij = Ao(7^j — 1) ~ Ao^i, and n = X^y/nijiji. The latter approx- 
imation is only chosen because it makes the linear response matrix appear simpler. The 
numerics in this paper are performed exactly, i.e., without this approximation. We divided 
the matrix into blocks. For very weak interaction strengths, the f-amplitudes are zero and 
the u-amplitudes are then solely determined by the upper left 2 x 2-block. The diagonals of 
this submatrix account for the excitation energy contributions due to the external and inter- 
action potential of the corresponding fragment. The off-diagonals account for the coupling 
energy to the other fragment. For stronger interaction strengths, the off-diagonal 2x2- 



blocks become important and induce finite w-amplitudes. As we have seen in Sec. IV, those 
lead for example to the damping of density oscillations. For spatially disjunct orbitals, i.e., 
j dx\(lPi^{x)\\(lPji{x) \ = 0, the linear response matrix boils down to two independent matrices, 
each acting on a separate subsystem. 
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We note that unlike the case M = 1, where the eigenvectors of 'P CfP can be obtained 
by apphcation of the projection operator on the eigenvectors of C [HO], this does not hold 
anymore for fragmented state^^ i.e., for M > 1. We demonstrate this property for M = 2. 
The eigenvectors of jC^^2 5 given by 



'M=2 







( VI) \ 






VI) 


\yt) 




K) 


V 1^^') ) 




\ \n) ) 



(B6) 



are in general not orthogonal to the ground-state orbitals. For the statement to be valid, 
the following expression must vanish: 



M=1 



{l-V) 



1 Vt) \ 


f 


V'r) 




\vt) 




\ K) 1 





mmu'R) + mmu'R) 



\\<pT){<i>T\yii) + \<PT){<i>T\yR)j 



0,*\ 



lO,* I 



(B7) 



This is not the case since the only eigenvector of 'PCj^^2 with eigenvalue zero which lies 
in the space spanned by the ground state vectors is 10?;), 10°'*), |0/j*))^j^ As stated 

above, the only exception is M = 1, where generally \u) = \PU) and \v) = \PV) hold. 



Then, by using the orthonormalization relations Eq. (37) for M = 1, we find that the 



vector appearing on the right hand side of Eq. (B7) turns out to be proportional to the 
zero-eigenvector (|0°), |0°'*)^) oi CbcLG [see Eq. Q] [HQ]. 



Appendix C: Comparison to the linear response of two-component GP equations 

In this Appendix we compare the linear response of a two-fold fragmented condensate 
derived here, see Appendix [b| to that of a two-component EEC, see for example Ref. [26] . 
The latter system is described by two coupled Gross-Pitaevskii (2GP) equations pH - [27] 

z0° (a;) = [hL{x) + AiiK - + ^LRnR\4>l{x)\^] 0° (a;) , 

= [hR{x) + \nR{nn - ml{x)\' + AiHnL|0i(x)p} 0?,(a;) . (CI) 

We note that it also does not hold for the multi-component GP equation, see Appendix jcj 

The other eigenvectors with eigenvalue zero can be constructed as vectors which are transformed by Cj^—^ 

into a linear combination of ground-state orbitals, similar as for BdG in Ref. |80) . 



37 



We borrow the "above" notation with 'L' and 'R' representing the two species. For each 
component, the single-particle Hamiltonian is given by /iL,i?(a^) = ~2rn^& ~^ Vl,r{x), 
taking into account the different mass and potential trap of each component. The equations 



are different from the TDMF equations [see Eq. (17)] in two respects. Firstly, the two- 



component GP equations (CI) do not contain projectors P, since the atoms in different 
components are distinguishable. Secondly, the factor of 2 which appears in the TDMF 
equations due to the exchange interactions between identical particles is absent here. Instead, 
for two component BECs we have three interaction parameters, where Xn and X^ir account 
for the interactions between atoms of the same species, and Xlr between atoms of different 
species. A dynamical comparison of single-component fragmented and two species BECs, 
examining the case of interaction-assisted self interference in free space, can be found in 
Ref. [59]. 

Those differences between the identical and distinguishable particles also translate to 
the linear response of the two-orbital TDMF and the 2GP equations. In particular, we 
have found that for fragmented single-species BECs the response is orthogonal to all of 
the ground-state orbitals. This is not, of course, the case with two coupled GP equations, 
where, even when one employs a number- conserving framework for linear response as in 
Ref. [SO], the excitations of a given species are orthogonal solely to the ground state of 
the same species. This becomes important, e.g., for two-component BECs in a double- well 
potential, where at sufficiently high barrier each component is localized in one well. Then 
the excitation of the left species in the right well can have the same (say Gaussian) shaped 
ground state as the other species. For identical bosons we found a different behavior, where 



all excitations had at least one node in order to ensure orthogonality, see Sec. IV Another 
property of the two-component case is that one can investigate observables depending only 
on one orbital (species) p5j, in contrast to identical atoms where this is not possible. 



Appendix D: Linear response of the Bose-Hubbard model 

The assumption underlying BMF is that the ground state of the double-well is a perfect 



mean-field fragmented state, i.e. a Fock state [see Eq. (51)]. We investigate here the effects 
of a small tunnel coupling due to overlapping BMF orbitals by employing the two-site Bose- 
Hubbard model [m [66] in order to describe the hopping between the two BMF orbitals 
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^-iid 4>r(^)- Hamiltonian is given by 

where and create an atom in the left and right locahzed orbitals, respectively. 
The tunnel coupling is given by Q(t) = 2 J dx(j)\{x)h{x,t)(j)\{x), the asymmetry by = 
/ dx(t)l{x)h{x)(t)l{x) - J dx(t)%{x)h{x)(t)%{x). k ^ dx\(j)l{x)\'^ ~ f / dx\(l)%{x)\'^ is the 
prefactor of the interaction term. H^^ acts on an (A^+l)-dimensional state vector C, which 
entries mark the probabilities for having n atoms in the left and — n in the right orbital. 

In order to get the hnear response of the BH model, which we call LR-BH, we employ 
a time-dependent perturbation of the external potential as 5h = f (x) sin (ut) . Since the 
resonance frequencies of the orbitals as obtained from LR-BMF are in general different 
than the energies for hopping, i.e., the resonance frequencies of LR-BH, we assume time- 
independent orbitals (p^{x) and (Pr{x). We note that a general study of the interplay between 
orbitals' and hopping excitations requires a full many-body analysis. 

The perturbation affects the tunnel couphng: Q{t) Q{t) + SQ{t) , with 

Sfl{t) ^2 sin {cut) Jdx(f>l{x)f{x)(f>°ji{x), (D2) 
as well as the energy difference between left and right states 

SE{t) = sin (cot) J dx {\<j>l{x)\' - \<Pl{x)\') fix) . (D3) 
Linearizing C C° -|- 5C, we obtain 

i5c = H^"sc + 6n{t) (^didji + a^a^) + 6E{t) (d{dL - aj^a^) C72 , (D4) 

which can be straightforwardly solved by 

5C = - sin (ut) J^^kC'/ioj -ojk) . (D5) 

k 

Hereby, C*^ and Uk correspond simply to the eigenstates and eigenfrequencies of H^^ , re- 
spectively. More importantly, the response weights 7fc = 7^ + Jk given by 

7fe" = {c'\dld^ + a;,ajc°)/2 j dx(t>l{x)f{x)4>l{x) , (D6) 



and 



Ik = (claia^ - at,a^|c°)/2 j dx {\ct>l{x)\' - \ct>l{x)\') f{x) . (d7) 
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For the response weights related to tunnehng, which are proportional to the overlap of the 
orbitals, we find for all setups discussed in this paper ^ 0. For the response weights 
related to the potential asymmetry, we observe that in a symmetric double-well tI' = for 
symmetric perturbations /(x). For asymmetric perturbations, 7I' is non-negligible only for a 
few low lying excitations of LR-BH, since it scales as ~ (C'^l |C°). For example, 

we checked that for = 100 only the lowest excitation of LR-BH gives rise to a nonzero 
response. Moreover, one can show that the total position-space density response [i.e., p{x,t) 
as in Eq. (45)] at the LR-BH resonance frequencies scales with ~ ^(C^|d|^a^ — a^a^|C°) 
and is thus marginal. 
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FIG. 5. (Color online) Position-space densities and density oscillations for the same double-well 
system as in Fig. [2| The interaction strength is Aq^V = 10 in the left panels, and XqN = 20 in 
the right panels. The ground state density of GP is shown in (a,b) by the broad gray solid line 
and is scaled for better comparison to the density of BMF plotted by the broad black dashed line. 
The GP and BMF results are seen to coincide. Gerade density oscillations of the indicated excited 
states are shown in (a-f). The LR-GP results are shown by the blue solid lines. The dashed red 
lines (orange dashed-dotted lines) show results of LR-BMF for direct (swapped) excitations. See 
text for more details. All quantities are dimensionless. 
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FIG. 6. (Color online) Momentum-space densities and density oscillations for the same double-well 
system as in Fig. [2] The interaction strength is AqA^ = 10. The ground state density of GP is shown 
in (a) by the broad gray solid line and is scaled for better comparison to the ground state density 
of BMF, which is shown by the broad black dashed line. In momentum space they are completely 
different. Momentum-space density oscillations of gerade type of the indicated excited states (see 
Fig. [5]) are shown in (a,d,g), and those of ungerade type in (b,e,h). Blue solid lines correspond 
to density oscillations of LR-GP and red dashed lines to direct excitations of LR-BMF. Swapped 
excitations of LR-BMF are shown in (c,f,i). Orange solid (dashed) lines correspond to gerade-type 
(ungerade-type) density oscillations. See text for more details. All quantities are dimensionless. 
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FIG. 7. (Color online) Example of position-space density oscillations for the same double-well 
system as in Fig.[2]and for XqN = 20. Gerade [see Fig. [5](f)] as well as ungerade density oscillations 
of LR-GP are shown by the blue solid and dashed-dotted lines, respectively. Direct gerade (swapped 
ungerade) density oscillations of LR-BMF are shown by the red dashed (orange dashed-double 
dotted) lines, and are multiplied by a factor of for easier comparison. See text for more details. 
All quantities are dimensionless. 
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FIG. 8. (Color online) Excitation spectra of a BEC in a symmetric double-well potential (a = 0) 
versus barrier height for a fixed interaction strength XqN = 20. The total atom number is = 100, 
and the overall harmonic confinement given by ujho = V2. We compare the linear response of 
LR-GP, shown by the blue solid lines, and LR-BMF, shown by the red dashed lines for direct, 
and orange dotted lines for swapped excitations. The excitations are grouped into pairs of lines 
with gerade-ungerade symmetry and marked with numbers. The swapped excitations of LR-BMF 
are marked with primed numbers. Inset: We plot the potential by the black solid line. The 
corresponding ground state orbital of GP for 6 = 14 is shown by the blue solid, and the left orbital 
of BMF by the red dashed line. All quantities are dimensionless. 
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FIG. 9. (Color online) Position-space densities and density oscillations for the same double-well 
system as in Fig. [8] for a low barrier height b = 14. The ground state density of GP is shown in 
(a) by the broad gray solid line and is scaled for better comparison to the BMF density shown 
by the broad black dashed line. The GP and BMF results are seen to almost coincide. The blue 
solid lines in (a,c,d) show density oscillations of LR-GP. The red dashed lines in (c,d) show density 
oscillations of direct excitations of LR-BMF. (a) Low-lying excitation of LR-GP, which is absent 
for LR-BMF. (b) The gerade and ungerade swapped excitations of LR-BMF, marked as 1' and 
shown by the orange solid and dashed lines, respectively, (c) The gerade and (d) the ungerade 
excitations of LR-GP and LR-BMF, marked as 1. See text for more details. All quantities are 
dimensionless. 
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FIG. 10. (Color online) Momentum-space densities and density oscillations for the same double- 
well system as in Fig. [8] for a low barrier height b = 14. The ground state density of GP is shown 
in (a) by the broad gray solid line and is scaled for better comparison to the one of BMF, which is 
shown by the broad black dashed line. We show the same indicated excited states as in Fig. |9j The 
blue solid lines in (a,c,d) show density oscillations of LR-GP. The red dashed lines in (c,d) show 
results for direct excitations of LR-BMF. (a) Low-lying excitation of LR-GP, which is absent for 
LR-BMF. (b) The gerade-type and ungerade-type swapped excitations of LR-BMF, marked as 1' 
and shown by the orange solid and dashed lines, respectively, (c) The gerade-type and (d) the 
ungerade-type excitations of LR-GP and LR-BMF, marked as 1. See text for more details. All 
quantities are dimensionless. 
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FIG. 11. (Color online) Excitation spectra of a BEC in an asymmetric double-well potential 
(a = 0.1) versus occupation of the left orbital. We choose a high barrier height 6 = 20 and 
interaction strength XqN = 10. The total atom number is = 100, and the overall harmonic 
confinement given by ooho = V^- The top of the filled gray area shows the BMF ground state 
energies as a function of occupation of the left orbital. The vertical dashed line marks the location 
of the optimal occupation ni/N 0.55. We compare the linear response of LR-GP, shown by 
the blue solid lines, and LR-BMF, shown by the red dashed lines for direct, and orange dotted 
lines for swapped excitations. The excitations are grouped into pairs of lines with gerade-ungerade 
symmetry and marked with numbers. The swapped excitations of LR-BMF are marked with 
primed numbers. Inset: We plot the potential by the black solid line. The corresponding ground 
state orbital of GP at the optimal occupation is shown by the blue solid. The left (right) orbital 
of BMF is shown by the red dashed (short-dashed) line. All quantities are dimensionless. 
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FIG. 12. (Color online) Position-space densities and density oscillations for the same double-well 
system as in Fig. [IT] at the optimal occupation ni/N 0.55. The ground state density of GP is 
shown in (a) by the broad gray solid line and is scaled for better comparison to the one of BMF, 
which is shown by the broad black dashed line in (b). The density oscillations for an indicated 



excited state show up two lines for LR-GP (see Fig. 11). The lower in energy density oscillations 

are shown by the blue solid, and those higher in energy by the blue dotted lines in panels (a,c,e). 

A similar statement holds for the direct excitations of LR-BMF, where the respective lower excited 

states are shown by the red solid, and the respective upper ones by the red dotted lines in panels 

(b,d,f). The lower (upper) swapped excitation are shown by the orange dashed (dashed-dotted) 

lines in (b,d,f). See text for more details. All quantities are dimensionless. 
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FIG. 13. (Color online) Momentum-space densities and density oscillations for the same double- 



well system as in Fig. 11 at the optimal occupation til/N ~ 0.55. The ground state density of 
GP is shown in (a) by the broad gray solid line and is scaled for better comparison to the one of 
BMF, which is shown by the broad black dashed line. We show the same indicated excitations as in 



Fig. 12 The density oscillations for a given number show up two lines for LR-GP (see Fig. 11 ). The 



lower in energy density oscillations are shown by the blue solid lines in (a,d,g), and those higher 
in energy by the blue dotted lines (b,e,h). A similar statement holds for the direct excitations of 
LR-BMF, where the respective lower excited states are shown by the red dashed lines in (a,d,g), 
and the respective upper ones in (b,e,h). The lower (upper) swapped excitations of LR-BMF are 
shown by the orange solid (dashed) lines in (c,f,i). See text for more details. All quantities are 
dimensionless. 
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